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We analyze the heating of interacting bosonic atoms in an optical lattice due to intensity fluctua- 
tions of the lasers forming the lattice. We focus in particular on fluctuations at low frequencies below 
the band gap frequency, such that the dynamics is restricted to the lowest band. We derive stochastic 
equations of motion, and analyze the effects on different many-body states, characterizing heating 
processes in both strongly and weakly interacting regimes. In the limit where the noise spectrum 
is flat at low frequencies, we can derive an effective master equation describing the dynamics. We 
compute heating rates and changes to characteristic correlation functions both in the perturbation 
theory limit and using a full time-dependent calculation of the stochastic many-body dynamics in 
one dimension based on time-dependent density-matrix-renormalization-group methods. 
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I. INTRODUCTION 

Ultracold atoms in optical lattices provide a clean and 
controllable realization of quantum dynamics of an iso- 
lated many-body system on a lattice [1-3]. The re- 
markable progress in optical lattice physics is underlined 
by recent experiments including the quantitative deter- 
mination of phase diagrams and critical phenomena of 
strongly interacting Hubbard models [4-11], studies of 
quantum magnetism [12] and nonequilibrium quench dy- 
namics [13-15]. A basic experimental challenge is the 
preparation of low entropy or low temperatures states 
in optical lattices, and to avoid possible heating mech- 
anisms [7, 13, 14]. Heating can either be due to fun- 
damental decoherence sources like spontaneous emission 
[16-19], or collisional losses [20-22], and also due to tech- 
nical noise, for example, amplitude or phase noise of the 
lasers [23] generating the optical lattice. While in a re- 
cent publication [24] we have described possible optical 
lattice schemes which are immune to laser intensity fluc- 
tuations, we will present below a detailed study of heat- 
ing of bosonic atoms in an optical lattice as a many-body 
non equilibrium problem. 

We will study below heating within a single hand Bose- 
Hubbard model where the tunneling and hopping param- 
eters are stochastic functions of time reflecting the in- 
tensity noise of the laser. We derive this model under 
the assumption that the noise spectrum contains only 
significant components below the band gap of the lat- 
tice, i.e. noise induced transitions to higher bands can 
be neglected. In our model the intensity fluctuations 
of the light act as a global noise^ which corresponds to 
the assumption that the spatial correlations of the laser 
fluctuations are certainly much larger than the size of 
the atomic cloud. The resulting stochastic schrodinger 
Equation for the Bose-Hubbard dynamics will be solved 
in detail in various limits and approximations. First, 
we will compute the heating rates and the time depen- 
dence of characteristic correlation functions in a pertur- 



bative calculation valid for short times. In the white 
noise limit for the intensity fluctuations we will be able 
to perform the stochastic average and derive a master 
equation for the many-particle systems containing both 
the Hubbard dynamics and the heating terms. In addi- 
tion, we will solve the stochastic many-body Schrodinger 
equation in a Gutzwiller mean field approximation, and 
in one dimension with a time-dependent density-matrix- 
renormalization-group (t-DMRG) technique as a multi- 
plicative stochastic differential equation. Besides com- 
puting the total average energy transfer to the system as 
part of the heating dynamics, we will also provide a de- 
tailed study of the excitations in the many body system 
as signatures of the applied noise. 

The paper is organized as follows. In Sec. II we de- 
rive a stochastic Schrodinger equation for cold atoms in 
an optical lattice in the presence of intensity noise. In 
Sec. II C we describe the methods we use to analyze the 
resulting nonequilibrium dynamics, including full time- 
dependent calculations based on t-DMRG calculations 
in one dimension (ID) and a Gutzwiller mean-field ap- 
proach for three-dimensional (3D) lattices. In Sec. HI 
we present the resulting time-dependent dynamics and 
discuss heating in different parameter regimes, and in 
Sec. IV we present a summary and outlook. 



II. STOCHASTIC MANY-BODY 
SCHRODINGER EQUATION 

In this section, we derive and discuss a stochastic 
many-body Schrodinger equation (SMBSE) for ultracold 
bosonic atoms in an optical lattice in the presence of in- 
tensity fluctuations of the laser generating the lattice. We 
are interested in a situation where the atoms in the opti- 
cal lattice are prepared in the lowest Bloch band with dy- 
namics described by a single band Bose-Hubbard model 
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FIG. 1: (Color online) (a) We consider bosons in an optical 
lattice where noise in the lattice depth leads to noise in the 
hopping amplitude and in the interaction energy, (b) Rel- 
ative change of the hopping and the interaction parameter 
with lattice depth in an isotropic three-dimensional cubic op- 
tical lattice, generated in the standard way by three counter- 
propagating laser beams. Note that the change in the hopping 
and interaction parameter is anticorrelated in this setup. 



tion of the above model starting from a multiband Hub- 
bard model, which establishes the validity of the SMBSE 
given above, and gives corrections due to interband tran- 
sitions. In Sec. II B we will discuss the derivation of 
a master equation for the averaged density operator 
= with ((...)) denoting a stochastic 

average over the noise. This is possible under the as- 
sumption of a white- noise limit, i.e., SV{t) is modeled by 
Gaussian white noise (within the single band model). In 
Sec. II G we will discuss mean- field and DMRG versions 
of the SMBSE, and their simulation. 

Finally we note that similar discussions can be found 
in work on lattice spectroscopy [26-31]. There the po- 
tential is not fluctuating stochastically but modulated 
periodically in time. 



[25], 

Here the hopping amplitude and on-site interaction en- 
ergy are denoted by J and U respectively. The oper- 
ators bi are the annihilation operators for particles at 
site i. For deep lattices simple arguments give the de- 
pendence J - ^ER{V/ERf^^exp{-2y^V/E^) and 

U - SEr{V/Er)^/^ on the depth V of the optical lat- 
tice (with Er being the lattice recoil energy). Thus an 
increase of the lattice depth V suppresses the tunneling, 
while at the same time the on-site interaction becomes 
larger (Fig. 1). The single-band tight-binding model is 
valid provided the on-site interaction and temperature 
are much lower than the gap to the first excited band. 

Laser intensity noise can be included in the Bose- 
Hubbard dynamics as V{t) = Vo^SV{t) with SV{t) fluc- 
tuations around the mean lattice depth Vq. For small 
fluctuations the tunneling J{t) = Jo + jySV{t) and on- 
site interaction U{t) = Uo-\-SV{t)^ will become stochas- 
tic variables (see Fig. 1), and we can write a stochastic 
many-body Schrodinger equation (SMBSE) {h = 1) 

> = -.(H(J„f/.) + H(^,f)m«))w 

= -i{H^H'SV{t))\^), (2.2) 

for a given noise model SV{t). The derivatives 
^ are evaluated at Vq such that H = H{Jo^Uo) and 
= H{jy,^) are time independent. While H in- 
duces coherent evolution according to the standard Bose- 
Hubbard Hamiltonian, describes the heating in the 
lowest lattice band due to intensity noise. We expect the 
above model to be valid provided the fluctuation spec- 
trum of SV{t) is narrow on the scale given the separation 
to the first Bloch band. Otherwise, the noise will excite 
atoms to the higher Bloch bands. 

While the above heuristic derivation is intuitively ob- 
vious, we summarize below in Sec. II A a rigorous deriva- 



A. Bose-Hubbard Model with Intensity Noise 

To derive Eq. (2.2) We consider atoms (of mass m) 
in an optical potential Vopt{x) = [Vq -\- 5V {t)] sin^ {kx) , 
which for simplicity of notation we assume to be one di- 
mensional. Here k is the wave vector of the laser generat- 
ing the lattice, which is related to the lattice constant a 
via a = Tv/k and sets an energy scale via the recoil energy 
Er = /c^/(2m). The full many body Hamiltonian, can 
be written in second quantization using the bosonic field 
operators iIj{x) ['^^(x)] that destroy (create) a particle at 
the position x as 

H{V{t)) = I dx^\x) [-^^^y(t)sm\kx)^ ^{x) 
+ dxilj^{x)ilj^{x)tlj{x)ip{x), (2.3) 

where g = ATih^as/m and is the 5- wave scattering 
length. 

We are interested in the limit where the stochastic 
5V{t) is much slower than the (fast) time scale asso- 
ciated with the band gap. For a given lattice depth 
V the Wannier states Wj^niXiV) form a complete ba- 
sis. Thus it is natural to employ an adiabatic (Born- 
Oppenheimer) picture, where field operators are ex- 
panded into instantaneous Wannier states, that is ijjix) = 
X^i,n ^3.n{x, V{t))hi^niy{t)), whcrc 6i,n(V'(t)) annihilates 
a boson in the instantaneous Wannier state Wj^ni^^ ^(^)) 
at site i in band n. In this way the single particle basis 
states keep track of variations of V{t) on a slow timescale. 
Non-adiabatic transitions due to the time dependence of 
the basis states are driven by fast changes in the lattice 
potential. As shown in Appendix A the coefficients of 
the wave function |^) in the time-dependent Fock ba- 
sis corresponding to this instantaneous Wannier states, 
({^i,n}|^)7 evolve according to 

4{Kn}i*) = {Mmv{t))+G{v{t))v{tm, 
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with 

n {ij) n,i 

+ ^EEf^{»}(^)^Ul^L^^»3&i,n4, (2.4) 

i {n} 
i,n-J,m 

(2.5) 

To simplify notation we suppressed the explicit depen- 
dence of the bosonic operators and the Fock basis states 
on the instantaneous lattices depth V{t). An indepen- 
dent derivation and discussion of this equation in the 
context of a deterministically modulated lattice depth 
can be found in [31]. 

Due to the localized nature of the Wannier func- 
tions only nearest-neighbor hopping and on-site inter- 
actions are considered in (2.4). The first term (2.4) is 
simply the multi-band Bose-Hubbard Hamiltonian with 
time-dependent parameters, corresponding to the time- 
dependent lattice depth. The second term in (2.5) arises 
from the time dependence of the basis we use to describe 
the system. It is proportional to the time derivative of 
the potential depth. In Appendix A we give a detailed 
discussion of this term and we show that it drives transi- 
tions between different bands, but does not couple states 
within the same band. More precisely, this interhand 
term mainly couples atoms to the second excited band 
such that the number of atoms in the lowest band Nq 
decreases in time as Nq/Nq ^ —27]^ 82- Here S2 is the 
noise spectrum at the transition frequency from the low- 
est to the second excited band and 77 = na^ja < 0.1 is 
the Lamb-Dicke parameter that compares the extension 
ao of the lowest band Wannier function to the lattice con- 
stant a. This sets the timescale on which the restriction 
to the lowest band is valid. For fluctuations that are slow 
on the time scale of the gap this term is off resonant and 
can be dropped. 

Under the same constraints also the first term in (2.4) 
can be restricted to the lowest band and by linearizing 
the dependence of J{V) and U (V) for small fluctuations 
on finds Eq. (2.2). In the following we will describe the 
heating dynamics on the basis of the single band model 
(2.2). 

B. White Noise Approximation and Master 
Equation 

We are interested in calculating the response of the 
many-body system to the noise stochastically averaged 
over the various realizations of SV{t). This stochastic 
averaging can be performed exactly if we make the as- 
sumption of white noise {{SV{t)SV{t'))) = SoS{t-t'). We 
note that this white noise approximation implies that the 
fluctuations of SV{t) are much faster than 1/J and 1//7, 



but are much slower than the transition frequency to the 
first excited band. We can thus write the SMBSE as a 
Stratonovich stochastic differential equation 

{S) d\^) = -iH\^)dt - iH'y^\^)dWt, (2.6) 

where 5*0 denotes the strength of the noise and dWt is 
a Wiener increment [32, 33]. This is a multiplicative dif- 
ferential equation, for which the averaging over the noise 
can be performed exactly to derive a (master) equation 

To derive the master equation for p, we find it conve- 
nient to transform the above equation to an Ito equation 

(/) d\^) = (^-iH - ^H'^^ \^!)dt - iH'^Q\^)dWt. 

Using the Ito rules for stochastic calculus [32, 33] one 
finds the following many-body master equation: 

^^p=-i[H,p\-^[H',[H',p\]. (2.7) 

We note that this equation is of Lindblad form. The first 
term on the right hand side is the familiar Bose-Hubbard 
Hamiltonian (2.1), while the second term describes heat- 
ing. We note that the assumption of global intensity 
noise is reflected in the spatially nonlocal heating terms 
contained in the double commutator. The above equa- 
tion is derived from averaging over classical noise (as op- 
posed to coupling to a quantum reservoir). As a con- 
sequence solutions p{t) will in general approach for long 
times p 1 corresponding to an (infinite temperature) 
completely mixed state (within the subspace allowed by 
the conserved quantities). 

In Sec. Ill below we will derive analytical, perturba- 
tive solutions of the master equation to describe the ini- 
tial heating of a many-body quantum state in the limit 
of weak (U <C J) and strong interactions {U ^ J) . 
However, for general parameters we find it more con- 
venient instead of solving the master equation numeri- 
cally to compute averages from simulating trajectories of 
SMBSE (c.f. Sec. HC). 

C. Simulation of the Stochastic Many-Body 
Schrodinger Equation 

A simulation of the SMBSE as a multiplicative stochas- 
tic differential equation can be performed in a mean-field 
limit and for ID systems using t-DMRG techniques [34- 
37] or exact state representation for small systems. 

The (mean field) Gutzwiller-ansatz [38-40] for the 
Bose-Hubbard model relies on a product state assump- 
tion 1^) = = III Z)n fi,n\n)i. The time-dependent 
variational ansatz [41, 42] for a homogeneous system 
leads to a nonlinear Stratonovich stochastic Schrodinger 
equation of the form 

{S) d\^i) = -iHi\^i)dt - i^QH[\c^i)dW, (2.8) 
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where 

Hi = -zJ (^f 6, + i,ib]) + ^bfbf, (2.9) 
HI = -z^ + i'M) + Y^bfbl (2.10) 

and = {(j)i\hi\(j)i) . The number of nearest neighbors is 
denoted by z. 

To simulate the stochastic differential equation we typ- 
ically used a semi-implicit method given in [33] of strong 
order 1.0. The evolution in a small time step At is cal- 
culated from 

|f = \^t) - '^H\^t)^t - '-^/SH'\^t)AW, 

\^t+At)^2\<S>t)-\^t), (2.11) 

with a randomly chosen Wiener increment chosen from a 
normal distribution AVI^ = W{At) - VI^(O) - A^CO, At). 
The density matrix and expectation values of operators 
are then obtained by averaging over the trajectories cal- 
culated in this way. We directly incorporate this prop- 
agation scheme for exact state representations and the 
analog version for the Gutz wilier equations (2.8). 

Similar techniques can be employed in t-DMRG. There 
it is more convenient to implement the propagation 
step in the form of a Trotter decomposition. There- 
fore, we write the ID Hamiltonian and noise term as 
sum over next-neighbor operators H = ^ • and 
= "^Zi^i respectively. For small time steps the 
evolution step can then be implemented as 

\%+At) « n e-^^*^'-+^ n e-^^^^-' (2.12) 

i i 

To lowest order this is equivalent to the Euler algorithm 
and of weak order 1.0 convergence [32, 33]. Note that we 
did not find any stability issues when using a number con- 
serving update of the matrix product states. We further 
confirmed that for sufficiently small timesteps the results 
from the propagation (2.12) coincide with the results ob- 
tained by exact state representation from (2.11) for small 
systems {N bosons on M sites with N = M < 10). 



III. NONEQUILIBRIUM MANY-BODY 
DYNAMICS AND HEATING 

In this section we show our results for the heating rates 
and analyze how the noise changes the characteristics of 
the many-body ground state in the system. We show how 
to obtain analytical results for heating rates in the two 
limiting cases of a weakly interacting condensate in the 
superfluid (SF) phase (/7 <C J) and of a nearly perfect 
Mott insulator (MI) {U ^ J) and compare them to nu- 
merical simulations in Sec. Ill A. In Sec. IIIB we analyze 
how the noise changes the characteristics of the state and 
analyze the evolution of the condensate fraction. 




FIG. 2: (Color online) Short-time heating rates of superfluid 
{U — 2 J) and Mott insulator states {U = 6 J) in one dimen- 
sion as a function of the relative magnitude of noise on J and 
U. We parametrize the correlations between the noise on J 
and U by and A as VS{dJ/dV)/J = -Acos^(6') for < 
< 7t/2] VS{dU/dV)/U = Xs'm^iO) and VS{dJ/dV)/J = 
Acos^(^) for 7r/2 < ^ < TT. The usual anticorrelated case cor- 
responds to < 7r/2, and the sweet spot ^ = of Eq. (3.2) 
to = 37r/4. The heating rates are calculated from linear 
regression over 500 t-DMRG trajectories in a system with 
30 particles on 30 sites (open boundary conditions). In both 
cases heating is strongly suppressed in the vicinity of the sweet 
spot. Thin lines show results in the presence of a harmonic 
trapwith£^/J = 0.0356^^75^^^ 5xlO"V"^/^; in both 

cases we used A = 0.02J~^/^, time step At = 10~^/J. 

A. Heating rates 

If the density operator is diagonal in the eigenstates 
of the Hamiltonian H = Hj + Hu (where Hj and Hu 
denote the kinetic- and interaction-energy terms in the 
Bose-Hubbard model), for example if the system is in 
the ground or a thermal state, the average increase of 
the energy E = (H) can be calculated from the master 
equation (2.7) and is 

The expectation value can be evaluated analytically in 
the limiting cases of an ideal superfluid state, and an 
ideal Mott insulating state, as discussed below. 
From (3.1) we see that the heating vanishes if 

ldJ_ldU 
JdV~UdV 

Only if this condition is met, the Hamiltonian H and the 
noise operator H' commute, moreover this means that 
they are proportional to each other. As a consequence 
all states that commute with such as energy eigen- 
states or thermal states, are stationary if (3.2) is satis- 
fied (see Fig. 2). In the standard setup, the hopping rate 
always decrease with the lattice depth, while the onsite 
interaction always increases, such that there is no such 
"sweet spot" [see Fig. 1(b)]. However, one can come up 
with more elaborate lattice setups (see, for example, [24]) 
that are designed in such a way to fulfill the sweet spot 



5 



(a) Eq/ER 




(b) 



FIG. 3: (color online) Elementary heating processes in the 
limiting cases of weak and strong interaction, (a) For weak 
interactions {U <^ J) lattice fluctuations create pairs of Bo- 
goliubov excitations with opposite quasimomenta on top of 
the Bose-Einstein condensate at g = 0. (b) In the strongly 
interacting case {U ^ J) lattice fluctuations create pairs of 
excess particles and holes of opposite quasimomenta zLq on 
top of the Mott insulator. 



condition and therefore are resilient against this type of 
noise. For later convenience we introduce the parameter 



that measures the deviation from this sweet spot. 



1 dJ 
JdV 



1. Weak interactions 

If the interactions are weak, the ground state of the 
system is a Bose-Einstein condensate, where a macro- 
scopic number of atoms occupies the mode with zero 
quasimomentum. For an ideal condensate we obtain from 
Eq. (3.1) (for a cubic lattice with z nearest neighbors and 
a filling of n atoms per site) to lowest order vnU/ J the 
heating rate per particle: 



m 

N 



JdV 



UdV 



zJU'^n. 



(3.4) 



To find corrections to this result we can apply a Bogoli- 
ubov approximation. This approximation is most conve- 
niently expressed in the Bloch basis rather than in the 
Wannier basis, such that the kinetic part of the Hamil- 
tonian reads Hj = ^qSqb^qbq, where the single-particle 
energy spectrum Eq for a cubic lattice in d dimensions 

is Eq = 2J^^^-l[1 — cos{qia)]. Here qi denotes the com- 
ponent of the quasi momentum along direction i. The 
interaction is treated on a mean field level, replacing 
Hu Un/2 T.q{2blbq^bqb-q^blblq). The total Hamil- 
tonian is then quadratic and can be diagonalized by a 
standard Bogoliubov transformation bq = UqCq + V-qC^_q 



such that H 



T.q£q4^q^ wlth 



= y/eq{eq^2Un). 
Within the same approximation the noise operator = 
jj^Hj -\- jj^Hu is also quadratic. However, the Bo- 
goliubov transformation diagonalizing H does not diag- 
onalize H' [except if condition (3.2) is met]. Therefore, 
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FIG. 4: (Color online) Comparison of the short-time heat- 
ing rates from the Gutzwiller ansatz and t-DMRG simula- 
tions to the analytical results for weak and strong interac- 
tions. The Gutzwiller calculation is for a homogeneous in- 
finite system, the DMRG simulation for a ID system of 30 
particles on 30 sites. We consider a small anticorrelated noise 
{^^j^ = -0.01J-^/^ ^^S^^ = O.OIJ-^/'). We aver- 
age over ut — 1000 {ut — 500 for t-DMRG) noise trajectories 
and estimate statistical errors of the mean. Convergence has 
been checked with time-steps 1 x 10~^ < AtJ < 1 x 10~^ 
and Hilbert space truncations 8 < < 16. The t-DMRG re- 
sults are converged with a bond dimension oi D — 200. The 
inset shows the averaged mean energy as a function of time 
(solid lines) together with the linear increase (dashed lines) 
with heating rate given by eq. (3.1). This shows that the en- 
ergy increase is well captured by a constant heating rate for 
several hopping times. The parameters and methods used to 
calculate the inset are the same as in Fig. 5(b). 

the noise operator contains terms cJcL^ and thus excites 
pairs of Bogoliubov excitations with opposite quasimo- 
menta [Fig. 3(a)]. The total quasimomentum is con- 
served as required by the conserved symmetry. The heat- 
ing rate associated with this process can easily be calcu- 
lated from (3.1) for the vacuum of quasiparticles. For a 
d-dimensional cubic lattice {z = 2d) it is given by 



N 



So 



1 dJ 
JdV 



1 dUY 
Udv) 



2JU^nG{nU/J), (3.5) 



with 



G(e) 



1 f , (Eti[l-cos(x,)]) 



3/2 



\<Xi<n (E^ll[l -C0S( 



1/2- 



(3.6) 



In one dimension this integral can be calculated exactly 
and leads to 



G{e) = 1 - ^ + -a/^ + -(2 - e) arcsm ---- 



(3.7) 



For arbitrary dimensions we can expand the integral and 
obtain (to first order in e) G{e) = d — |, such that for 
U/J^O the rates (3.4) and (3.5) coincide. 
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2. Strong interactions 



3. Numerical results 



According to Eq. (3.1), the heating rate for a perfect 
Mott insulator with an integer number of atoms per site 
n is (to lowest order in J /U) given by 



N 



Also here we can obtain corrections as well as insight 
into the excitation process by suitable approximations. 
Following [44] one can approximately describe the one- 
dimensional Bose Hubbard model in the Mott regime by 
restricting the local Hilbert space to the three states \n) 
and |n±l). Using a generalized Jordan-Wigner transfor- 
mation one introduces fermionic creation operators for 
the excess particles (c]+) and holes (cj_). Assuming 
that the density of excess particles and holes is small, 
the Hamiltonian H = Hj + Hu and the noise operator 
H' = j^Hj + jj^Hu can be written in the quasimo- 
mentum basis approximately as 

Hj ^ -2J^cos(pa)[(n + l)cj^^^Cp,+ nc\^_Cp-] 



+ V'n(n + l)zsin(pa)(cj^^^cLp^_ 
Hu ^ncji^+Cp,+ - (n - l)cl^_ 

k 



-), 



(3.9) 



The Hamiltonian H is then quadratic and can be di- 
agonalized by the Bogoliubov transformation = 

UpCl^^ + VpC^p^^, such that H = Z)p,a=± ^^c^iphl^alp^cr 

with the quasiparticle dispersion relation 
e±(p) = -Jcos{pa) + -^(271 - 1) 



±-\J[U - 2(2n + 1) Jcos(pa)]2 + 16n(n + 1) sin^(pa). 

(3.10) 

As in the weakly interacting case, the noise operator 
is in general not diagonal in the quasiparticle basis, but it 
contains terms jp^_^j^_p^_ that generate correlated quasi- 
particle pairs that travel though the system with opposite 
quasimomentum. To lowest order in J//7, these are pairs 
of excess atoms and holes [Fig. 3(b)]. The heating rate 
associated with these processes can easily be calculated 
from (3.1) with (3.9) and (3.10) for the vacuum of quasi- 
particles. It is given by 

2 



N 



So 



1 dJ 
JdV 



1 dU 
UdV 



where we abbreviated 



F(j,n) 



1 



dp 



2 sin^ (p) 



2Ur (n + l)F(J//7,n), 
(3.11) 



(3.12) 



/(j, n,p) = (1 — (4n + 2)j cos(p)) + 16(n^ + n)j^ sin^(p). 

For j <Cl one finds F(j,n) = l^{l-2n-2n'^)f ^0{j^). 
For J/U this reduces to (3.8). 



In Fig. 4 we show results for the short time mean heat- 
ing rate per particle for unit filling (n = 1) at differ- 
ent values for the interactions. The analytical results for 
weak and strong interaction are shown together with the 
results from numerical methods outlined in Sec. II C. The 
heating rates obtained with t-DMRG in one dimension 
(up to N = M = 30 particles per lattice site) agree very 
well with the analytical results in both limiting cases and 
connect these smoothly across the phase transition. 

On a mean-field level the (ground state) phase tran- 
sition from a SF (IV^^P > 0) to an MI {ipi = 0) ground 
state occurs at Uc ~ b.8z. The Gutzwiller wavefunction 
captures the two limiting cases of an ideal superfiuid as 
a product of coherent states at each lattice site for van- 
ishing interaction and a Mott insulator as a product of 
Fock states at each lattice site. However, a general lim- 
itation of the Gutzwiller mean-field theory is that the 
entire Mott insulating phase (at integer filling n) is rep- 
resented by the same wavefunction |^) = Yii This 
state is trivially invariant under the evolution with the 
stochastic equations (2.8). As an (unphysical) artifact of 
this limitation Eqs. (2.8) predict no heating in the entire 
Mott phase. The only nontrivial dynamics can be ob- 
served on the superfiuid side of the phase transition. On 
this side, except close to the phase transition, where the 
mean-field treatment is expected to fail, the results are in 
very good agreement with the analytical results obtained 
from the Bogoliubov approximation. 



B. Evolution of state characteristics 

The heating of a general many-body quantum state 
cannot be fully understood by a single heating rate, since 
the system is driven out of thermal equilibrium in gen- 
eral. In order to further quantify the heating, we analyze 
how characteristic correlation functions of the different 
many-body states are affected by the noisy lattice. 



1. Single Particle Density Matrix and Condensate fraction 

For the Bose-Hubbard model, the MI and the SF states 
are characterized by the off-diagonal correlations, i.e., the 
off-diagonal elements of the single-particle density matrix 
(SPDM), (blbl^j). The signature of the SF ground state 
is off-diagonal long-range order, i.e., these elements de- 
cay to a constant (decay algebraically in one dimension) , 
whereas they decay exponentially to zero in the MI. Here 
we analyze how these characteristics change as a function 
of time. 

Closely related to the SPDM is the condensate frac- 
tion, which is defined as the largest of the eigenvalues 
{AJ of the SPDM: J=' = Xo/J^i^i- ^o^e that in the 
case of the Gutzwiller ansatz, the condensate fraction 
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FIG. 5: (Color online) (a) The time evolution of the condensate fraction for short times in a large system with M = A/" = 48 
(averaged over 60 trajectories, bond dimension D = 256, local dimension truncation di = 6), (b) for long times in a small 
system with N = M = 10 (500 trajectories), and (c) for a homogeneous infinite Bose-Hubbard model in mean- field theory 
(Gutzwiller ansatz, di = 8, 1000 trajectories, solid grey lines obtained by a small noise approximation [33]). The fraction in 
general decreases except for the extreme Mott insulating case oi U / J — 30 and in the long time limit in mean-field theory. 
In (a) the noise is anticorrelated with v^jl^ = -0.025 J"^/^ = 0.025 J"^/^ in (c) with = -0.01J"^/^ 

^/Sjj^ — O.OIJ"^^^. In all three plots, the straight dashed green lines correspond to the result obtained in perturbation 
theory (see text). 



is simply given hy F = = \{^i\hi\(l)i)\^ . The per- 

turbative analysis Sec. Ill A shows that the main effect 
of noise in the lattice depth on a Bose-Einstein conden- 
sate is the generation of quasiparticle pairs, and there- 
fore a decrease of the number of atoms in the conden- 
sate mode with time. As for the heating rate one can 
calculate this depletion rate from the master equation 

(2.7) . Denoting the number of atoms in the mode with 
quasi-momentum q hy Nq = {^\^q) we have {{Nq)) = 
^^Uwd i^u^ blbq] , Hu] ) for the system initially in the 
ground state. In the limit of an ideal condensate the num- 
ber of atoms in the condensate mode ai q = evaluates to 
{{No))/No = -So^^^W^^. Within the same approx- 
imation the quasimomentum of the particles scattered 
out of the condensate mode is distributed homogeneously 

over the whole Brillouin zone: {{Nq)) = Sq(^^%U^^. 

Numerical solutions of the stochastic differential equa- 
tion show such a depletion of the condensate mode, both 
using exact methods [Figs. 5(a), 5(b) and Fig. 6] (exact 
diagonalization and t-DMRG) as well as in the Gutzwiller 
framework [Fig. 5(c)]. By performing a small noise ex- 
pansion [33] of the nonlinear stochastic Gutzwiller equa- 
tions we find that this initial decay of the condensate 
fraction in the Gutzwiller framework is associated with 
the incoherent excitation of the amplitude mode at zero 
quasimomentum [43] [see Fig. 5(c)]. However, we note 
that for long times solutions to the Gutzwiller equations 

(2.8) have a mean condensate fraction of approximately 
~ 0.3. This steady-state mean condensate fraction is 
not associated with a steady-state condensate fraction in 
the individual trajectory, where in the long-time limit 
the condensate fraction oscillates between and 1 with 
a random phase, which gives rise to a nonzero steady- 
state behavior for the stochastic average. Such a signa- 
ture is considered an unphysical artifact of the incapacity 
of the Gutzwiller wavefunction to capture the decay of 
the amplitude mode. In fact, exact simulations of one- 
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FIG. 6: (Colour online) The evolution of the population in 
different quasi momentum modes. Shown are averages over 
500 noise realizations for a superfluid state {U I J — 2) for 
a ID system with N — \0 particles on M = 10 sites. The 
lattice fluctuations lead to a decrease of the population in the 
condensate mode at g = in agreement with the prediction 
from perturbation theory given by the dashed line. The noise 
is anticorrelated with V^jJ^ = -0.025 J"^/^ = 
0.025J"^/^ 



dimensional systems show no such behavior. 

In Fig. 5(b), we show the stochastic average of the con- 
densate fraction for a one dimensional system of ten sites 
and periodic boundary conditions using exact diagonal- 
ization. Even though in one dimension the ground state 
shows only quasi- long-range order, for small values of 
U I J the decay of the condensate fraction is well described 
by the expression obtained in perturbation theory. Our 
t-DMRG simulations essentially lead to the same results 
but are limited to very short times [Fig. 5(a)]. In con- 
trast to the Gutzwiller results, the condensate fraction 
decreases monotonically also for long times for all super- 
fluid initial states. 
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FIG. 7: (Color online) (a) Parity correlations in the ground 
state of a small system with N — M — 10. The next- 
neighbor and longer-ranged correlations assume a maximum 
at t// J ~ 6 (solid lines are for open, dashed lines are for peri- 
odic boundary conditions), (b) Long-time evolution of these 
correlations calculated with exact diagonalization in the small 
system, averaged over 500 noise trajectories (periodic bound- 
ary conditions). The transient state develops long-range cor- 
relations, (c) The evolution of the parity-parity correlation 
function of a single noise trajectory in a Mott insulating state 
with U/J = 6 in a ID system with AT = M = 48 (bond di- 
mension D — 265, local dimension truncation di =6). The 
amplitude noise excites single particle-hole pairs which spread 
out through the system as a light cone. In all simulations, 
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the noise is anticorrelated with j ^ 

^ ^ U dV 



2. Particle-hole correlations 

In the limit of strong interactions, as shown above, 
we expect elementary excitations to consist of correlated 
particle-hole pairs propagating through the system. To 
analyze the dynamics associated with these excitations, 
we compute parity correlation functions defined as [44] 



(3.13) 



where sj is the local parity operator at site j, defined 
as Sj = exp{z7r(nj — n)}. Since we calculate these func- 
tions in finite inhomogeneous ID systems, depends 
on the site i from which we start calculating the function, 
and we will typically begin from the central site i = M/2 
in a system of M sites. In Fig. 7 we plot the evolu- 
tion of this correlation function under the evolution with 
the SMBSE, calculated with t-DMRG for a large system 
and with exact diagonalization for a small system. Note 
that useful information about these correlations cannot 
be obtained from a Gutzwiller product state ansatz, as 
the first term of (3.13) will always factorize and thus the 
correlation will always be zero. 

As seen in Fig. 7(a), the parity correlations in the 
ground state assume the largest value in the Mott in- 
sulating phase at /7/J ~ 6. In a noisy time evolution 



we find that the initially large next-neighbor parity cor- 
relation starts to decrease, whereas the long-range parity 
correlations start to increase, as shown in Fig. 7(b). It 
reaches a maximum at a transient state at times t J ~ 30, 
where the next-neighbor and long-range interactions as- 
sume nearly the same value. At longer times all these 
correlations start to decrease into a more and more clas- 
sical state without correlations. Furthermore, looking 
at a single noise trajectory in a large system, we find 
in Fig. 7(c) that the elementary excitations induced by 
the amplitude noise are seen as excitations in the parity- 
parity correlations, which spread out in the form of a 
light cone, similar to the results in [44], where these cor- 
relation functions are directly measured in a quantum 
gas microscope experiment. 



IV. SUMMARY AND OUTLOOK 

In conclusion we have derived a microscopic model for 
the dynamics of bosonic atoms in a time-dependent op- 
tical lattice. While controlled periodic modulation of 
the lattice has the potential to access interesting physics 
[29, 30], we considered the alternative situation where the 
lattice depth fluctuates stochastically at low frequency. 
This situation arises naturally in optical lattice experi- 
ments due to intensity fluctuations of the lattice lasers. 
Using analytical approximations as well as t-DMRG and 
time-dependent Gutzwiller methods, we have analyzed 
the nonequilibrium dynamics of many-body states in a 
variety of parameter regimes. We find characteristic re- 
sponses of the system that vary in different parameter 
regimes, and could be used to identify the effects of such 
dynamics in current experiments. 

The generalization of the initial model to fermionic 
atoms results in a stochastic equation of motion similar 
to (2.6), 

{S) d\^) = -iH\^)dt - iH'y^\^)dWt, (4.1) 

with a coherent part given by the Fermi-Hubbard Hamil- 
tonian, and a corresponding stochastic contribution: 

H = -J + ^ E 4t'^'.tcUcu, (4-2) 

^' = -^ E 4<.ci.. + ^E4tc^.t4iC^4- (4-3) 



Here the Ci^a are operators annihilating a fermion at site 
i with spin a G {t^i}- An analysis of this equation, 
similar to the one presented in Sec. Ill, can be carried 
out. In the Mott insulator at half filling for the anti- 
ferromagnetic ground state (/7 ^ J) a heating rate of 
E/N = Sq^^^UzJ'^ can be found from mean-field calcu- 
lations. Exploration of heating in other regimes, e.g., in 
the BCS-BEC crossover, is an interesting direction for 
further analysis. 
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term couples only states in different bands it is instruc- 
tive to transform G{V) to the instantaneous Bloch basis 
4^q,n{x^ y) (with the corresponding annihilation operator 

(A4) 

The derivative of the Bloch state (j)q^ri{x^ V) with respect 
to the lattice depth can be expressed in terms of the Bloch 
states with the same quasimomentum but in different 
bands m ^ n. From ffist-order perturbation theory we 
find 



Appendix A: Interband transitions 

Here we comment on the derivation of the SMBSE fo- 
cusing on the transitions to higher bands due to nonadi- 
abatic transitions. 

To proceed from Eq. (2.3) we expand the field operator 
in the instantaneous Wannier basis. The coefficients of 
the wavefunction |^) in the corresponding instantaneous 
Fock basis |{^i,n}) change in time via the change of the 
wavefunction and via the change of the time-dependent 
basis states: 



dt 



(Kn}|^) = 



|(K.}|)w + (K4|(||^) 



(Al) 



The time derivative of a Fock basis state |{^i,n}) 
constructed from the time dependent single particle 
basis Wi^ri{x^V{t)) can be obtained from the time 
derivative of the corresponding annihilation opera- 
tors: bi^n = J dxWi^n{x^V{t))ilj{x). With ?/^(x) = 
5^2,71 ^^,n(^,V'(t))6i,n we find 



d 

dt 



dx ( ^^Wi^n{x,V{t)) ) Wi^n{x,V{t))hi 



(A2) 



This expresses the change of the annihilation operators 
in the Schrodinger picture due to the change of the basis 
states. It should not be confused with the Heisenberg 
equation of motion for this operators. With this relation 
it is easy to show that the corresponding Fock states 
|{^i,n}) change due to the changing single particle basis 
as 



d 
dt 



IK,n}) 



J dxWj^rn{x,V{t)) (^-^Wr,s{x,V{t))^ ^j,m^r,s | {^i,n}) 

(A3) 



which gives rise to the term G{V{t)) in Eq. (2.5). A 
similar derivation is given in [31]. To see that this 



dV 



,,n(^,V) = ^ (l)q,m{x,V)aq^rn{V), 



^g,m(V) 



/ dycl^lmiy. V)0g,n(^, V) sin^(%) 



(A5) 



(A6) 



Thus the integral in (A4) is nonzero only for p = q and 
for n 7^ m. The operator G{V) therefore transfers sin- 
gle atoms into a different band with the same quasimo- 
mentum. This conservation of quasimomentum reflects 
the fact that the lattice translation symmetry is not bro- 
ken by global fluctuations in the lattice depth. Further, 
the perturbation does not break the reflection symmetry 
X —X. As a consequence, bands with even (odd) band 
index n are coupled only to bands with an even (odd) 
band index m. 

To lowest order in J/\oJn,m\^ with uJn,m being the en- 
ergy difference between the bands n and m, we can find 
a simple expression for the resonant parts in G{Vo)SV{t) 
that affect atoms in the lowest band: 

G{Vo)SV{t) « Yl KnSV^{t)blji,o + H.c, (A7) 

where we again kept only terms diagonal in the site in- 
dex i due to the localized form of the Wannier functions. 
Here 5Vn{t) denotes restriction of the noise term SV{t) to 
frequencies around transition frequency from the lowest 
to the nth band. Approximating the Wannier functions 
at a given lattice site with harmonic oscillator wavefunc- 
tions, we find that atoms in the lowest band are dom- 
inantly scattered into the second excited band in such 
band-changing processes. To lowest order in the Lamb 
Dicke parameter rj = /cao, which measures the extension 
of the Wannier function ao, compared to the lattice con- 
stant a = 7r//c, one finds Ai,n ^ V2r]'^Sn,2^0{r]^). Thus, 
the number of atoms in the lowest band A^o decreases in 
time as Nq/Nq = - En>o ^nA,?^ ^ -2r]^S2- This sets 
the timescale on which the restriction to the lowest band 
is valid. 
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